probit.m
This script calculates population values of Bresnahan and Reiss's estimates of
based on their ordered probit model of entry.
Contents
- Undertake the prerequisite calculations, if necessary.
- Create an in-line function for the standard normal c.d.f.
- Eliminate all sates with near-zero probability.
- Set up the binary likelihood function and its derivatives.
- Maximize the binary likelihood functions.
- Set up the full multinomial likelihood function.
- Maximize the ordered probit likelihood function.
- Report results if this is a test run.
Undertake the prerequisite calculations, if necessary.
The required inputs are CN, the matrix of possible (C,N) pairs and v, the vector with each pair's probability in the ergodic distribution.
if ~exist('probitARG','var') bellman thresholds ergodic CN=[kron(ones(maxN-minN+1,1),omega') kron((minN:1:maxN)',ones(nCstates,1))]; else CN=probitARG.cn; v=probitARG.v; end CData = exp(CN(:,1)); NData = CN(:,2);
Current plot held Current plot held Current plot held Current plot held

Create an in-line function for the standard normal c.d.f.
Phi = @(x) 0.5*(x>=0).*(erf(abs(x)/sqrt(2))+1)+0.5*(x<0).*erfc(abs(x)/sqrt(2));
Eliminate all sates with near-zero probability.
This avoids numerical errors arising from multiplying very, very small probabilities by ``infinite'' log likelihood evaluations.
eps=1e-10; CData=CData(v>eps); NData=NData(v>eps); v=v(v>eps); v=v/sum(v);
Set up the binary likelihood function and its derivatives.
We first estimate the entry thresholds with standard binary probits. These estimates then serve as starting values for the ordered probit estimation. So that it can be used for all of the thresholds, the rank under consideration is read as one of the likelihood function's inputs.
lf = @(x) v'*((NData>=x(3)).*log(Phi(log(CData/x(1))/x(2))) + (NData<x(3)).*log(Phi(log(x(1)./CData)/x(2))));
% Calculate the minimum and maximum values of N consistent with the ergodic distribution.
Nlow=min(NData);
Nhigh=max(NData);
Maximize the binary likelihood functions.
This uses the optimization toolbox for its minimization.
bProbitC=zeros(Nhigh-Nlow-1,1); bProbitSigma=bProbitC; for R=Nlow+1:Nhigh; f = @(x) -lf([x R]); x=fminsearch(f,[0.5*(exp(overlineC(R))+exp(underlineC(R))) 0.2]); bProbitC(R)=x(1); bProbitSigma(R)=x(2); end
Set up the full multinomial likelihood function.
As written, x must be a row vector with [sigma epsilon Cbar(Nlow) Cbar(Nlow+1) ... Cbar(Nhigh) infinity].
logprobs = @(x) log(Phi(log(CData./x(NData-Nlow+2)')/x(1))... -Phi(log(CData./x(NData-Nlow+3)')/x(1))); %Likelihood function. First argument is standard deviation, and remaining arguments are thresholds. (``zero'' and %``infinity'' coming first and last). lf = @(x) v'*logprobs([ x(1) 1e-10 x(2:end) 1e10]);
Maximize the ordered probit likelihood function.
f= @(x) -lf(x); xm=fminsearch(f,[mean(bProbitSigma) bProbitC']); mProbitC=xm(2:end); mProbitSigma=xm(1);
Exiting: Maximum number of function evaluations has been exceeded - increase MaxFunEvals option. Current function value: Inf
Report results if this is a test run.
if ~exist('probitARG','var') disp(mProbitC); disp(mProbitSigma); disp(f(xm)); end
0 0.7164 1.1383 2.0959 0.3631 Inf